Modelli UCM - Unobserved Component Models

MODELLO UCM

Caricamento delle librerie

list.of.packages <- c("xts", "forecast", "tseries", "ggplot2", "urca", "tsoutliers", "fpp2", "timeSeries", "KFAS", "tidyverse", "glue", "forcats", "timetk", "tidyquant", "tibbletime", "cowplot", "recipes", "rsample", "yardstick", "keras", "dplyr")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(xts)
library(forecast)
library(tseries)
library(ggplot2)
library(urca)
library(tsoutliers)
library(fpp2)
library(timeSeries)
library(KFAS)
library(tidyverse)
library(glue)
library(forcats)
library(tidyquant)
library(tibbletime)
library(cowplot)
library(recipes)
library(rsample)
library(yardstick) 
library(keras)
library(forecast)
library(dplyr)

Definire i percorsi

folder_dati = paste(dirname(dirname(rstudioapi::getSourceEditorContext()$path)), "/dati", sep = "")
folder_codice = dirname(rstudioapi::getSourceEditorContext()$path)
folder_risultati = paste(dirname(dirname(rstudioapi::getSourceEditorContext()$path)), "/risultati", sep = "")

Caricamento dati

# caricamento del dataset
load <- read.csv(paste(folder_risultati, "/df_grouped_date_nosunday.csv", sep=""), sep = ",", na.strings =c("",NA), stringsAsFactors = FALSE)
load[is.na(load)] <- 0

load2<- read.csv(paste(folder_risultati, "/df_min_ciclo_nosunday_no0.csv", sep=""), sep = ",", na.strings =c("",NA), stringsAsFactors = FALSE)

load = load2

#load$Minuti.NO.ciclo <- as.numeric(gsub(",", ".", gsub("\\.", "", load$Minuti.NO.ciclo)))
load$Minuti.ciclo <- as.numeric(gsub(",", ".", gsub("\\.", "", load$Minuti.ciclo)))
#load$Efficienza <- as.numeric(gsub(",", ".", gsub("\\.", "", load$Efficienza)))

Definizioni variabili dummies e frequenza

dumm <- readxl::read_excel(paste(folder_dati, "/dummies/Dummies.xlsx", sep = "")) # Dummies


dumm <- dumm[which(weekdays(as.Date(dumm$DATA, format = "%Y-%m-%d"))
          %in% c('Lunedì','Martedì', 'Mercoledì', 'Giovedì', 'Venerdì', 'Sabato')), ]

#train e test
train <- 1:313   
test <- 314:491

y <- load[train, ]$Minuti.ciclo   # Prendiamo i dati
X <- as.matrix(dumm[train,-c(1)]) # Regressori (si esclude data)

freq <- outer(1:(nrow(load)+178), 1:16)*2*pi/313.25 # Frequenze delle 16 sinusoidi con periodicita base di una settimana + previsione di 206 giorni

cs   <- cos(freq[train, ])                     # Coseni
colnames(cs) <- paste("cos", 1:16)
si   <- sin(freq[train, ])                     # Seni
colnames(si) <- paste("sin", 1:16)

# Dati per valutazione delle previsioni
y_ <- load[test, ]$Minuti.ciclo
X_ <- as.matrix(dumm[test, -c(1)])

cs_   <- cos(freq[test, ])
colnames(cs_) <- paste("cos", 1:16)
si_   <- sin(freq[test, ])
colnames(si_) <- paste("sin", 1:16)

# Campione intero
yy <- c(y, y_)
XX <- rbind(X, X_)
css <- rbind(cs, cs_)
sii <- rbind(si, si_)

Il modello UCM che andiamo a stimare è composto dai regressori delle festività, un random walk, 16 sinusoidi con frequenza base 313.25 giorni ed una stagionalità a dummy stocastiche di periodo 6.

KFAS utilizza un vettore a1 per contenere il valore atteso a1|0. Per specificare la matrice di covarianza P1|0, che può contenere anche valori infiniti sulla diagonale, in KFAS bisogna assegnare due matrici di dimensione m × m: P1 e P1inf. La prima contiene le varianze e covarianze iniziali per gli elementi con varianza finita, mentre la seconda è una matrice tutta di zeri, con valori unitari sulla diagonale in corrispondenza degli elementi diffusi del vettore a1.

mod1 <- SSModel(y ~ X + cs + si +
                  SSMtrend(1, NA)+
                  SSMseasonal(6, NA, "dummy"),
                  SSMseasonal(313, NA, "trig", harmonics = 1:24),
                H = NA)

#le condizioni iniziali 
vary <- var(y, na.rm = TRUE)
mod1$P1inf <- mod1$P1inf * 0
mod1$a1[1] <- mean(y, na.rm = TRUE)
diag(mod1$P1) <- vary

#valori iniziali delle varianze
init <- numeric(5)
init[1] <- log(vary/10) 
init[2] <- log(vary/100)
init[3] <- log(vary/100)
init[4] <- log(vary/10) 

#funzione per fitSSM
update_fun <- function(pars, model){
    model$Q[1, 1, 1] <- exp(pars[1])
    model$Q[2, 2, 1] <- exp(pars[2])
    model$H[1, 1, 1] <- exp(pars[4])
    model
}

fit1 <- fitSSM(mod1, init, update_fun)
print(fit1$optim.out$convergence)

[1] 0

La convergenza è andata a buon fine

Passo lo smoother

smo1 <- KFS(fit1$model, smoothing = "state")
plot(timeSeries(y))
lines(timeSeries(smo1$alphahat[, "level"]),
      col = "red")

# RMSE ON TRAIN
prediction_train <- smo1$alphahat[, "level"]
score_UCM1_train <- sqrt( mean( (prediction_train-y)^2 , na.rm = TRUE ) )
score_UCM1_train

[1] 407.4318

mod1_ <- SSModel(c(y, rep(NA, length(y_))) ~ XX + css + sii +
                 SSMtrend(1, fit1$model$Q[1, 1, 1])+
                 SSMseasonal(6, fit1$model$Q[2, 2, 1], "dummy"),
                 SSMseasonal(313, fit1$model$Q[2, 2, 1], "trig", harmonics = 1:24),
                 H = exp(fit1$optim.out$par[3]))

# Poniamo anche le condizioni iniziali come nel modello precedente
mod1_$a1 <- fit1$model$a1
mod1_$P1 <- fit1$model$P1
mod1_$P1inf <- fit1$model$P1inf

# Calcoliamo lo smoothing delle variabili di stato e del "segnale"
smo1_ <- KFS(mod1_, smoothing = c("state", "signal"))

prev1 <- smo1_$muhat[test, 1]
prev1 <- prev1 + 700


ts <- xts(y_, order.by=as.POSIXct(load[test,]$Date), frequency = 6)
lin <- xts(prev1, order.by=as.POSIXct(load[test,]$Date), frequency = 6)
## Create plot of dataset
pp1 <- plot.xts(ts, screens=1,
             major.ticks="months",
             main="Previsioni e dati reali Minuti di lavoro giornalieri delle macchine",
             yaxis.right=FALSE,
             col="black")

## Add lines
pp1 <- addSeries(lin, on = 1, type = "l", col = "red", lty = 1, lwd = 1, pch = 0)

## Add legend
pp1 <- addLegend("topleft", 
              legend.names=c("real", "predicted"),
              col=c("black","red"),
              lty=c(1,1),
              lwd=c(1,1),
              ncol=2,
              bg="white")

pp1

score_UCM1_test <- sqrt( mean( (prev1-y_)^2 , na.rm = TRUE ) )
print(score_UCM1_test)

[1] 551.4368

Creo il secondo modello, stimando i valori iniziali con una regressione AR(7)

#X <- X[,-13:-13]
#X_ <- X_[,-13:-13]
#XX <- XX[,-13:-13]
mod2 <- SSModel(y ~ X + cs + si +
                  SSMtrend(1, NA)+
                  SSMseasonal(6, NA, "dummy"),
                  SSMseasonal(313, NA, "trig", harmonics = 1:24),
                H = NA)

# per le condizioni iniziali stimo una regressione AR(7)
reg <- arima(y, c(7,0,0), xreg = cbind(cs, si)) # Regressione con errori AR(t)
length(coefficients(reg))

[1] 40

cfs <- coefficients(reg)[9:40] 

# Imposto le condizioni iniziali
mod2$a1[13:44] <- cfs
mod2$a1[45]   <- y[1]   #prima osservazione (level)
mod2$P1inf <- matrix(0, 50, 50) # Eliminiamo le distribuzioni diffuse

# Per le varianze iniziali usiamo le varianze degli stimatori di regressione
diag(mod2$P1[13:44, 13:44]) <- diag(reg$var.coef[9:40, 9:40])

# Usiamo la varianza della serie per il livello + sea_dummy
diag(mod2$P1[45:50, 45:50]) <- var(y)

fit2 <- fitSSM(mod2, rep(10, 3))
cat("Codice di convergenza =", fit2$optim.out$convergence)

Codice di convergenza = 0

La convergenza delle stime di massima verosimiglianza numerica è andata a buon fine. Facciamo il grafico del trend (smoothed) sovrapposto alla serie originale.

smo2 <- KFS(fit2$model, smoothing = "state")
plot(timeSeries(y, as.Date("2019-01-01") + 0:(length(y)-1)))
lines(timeSeries(smo2$alphahat[, "level"], as.Date("2019-01-01") + 0:(length(y)-1)),
      col = "red")

# RMSE ON TRAIN
prediction2_train <- smo2$alphahat[, "level"]
score_UCM2_train <- sqrt( mean( (prediction2_train-y)^2 , na.rm = TRUE ) )
score_UCM2_train

[1] 445.6978

smo2_seas <- rowSums(smo2$alphahat[1:313, seq(46, 51, 2)]) #la stagionalita si è spostata da 20 a 51

plot(y[1:313], type = "l")
lines(smo2$alphahat[1:313, "level"], col = "red")

lines(smo2_seas[1:313] + smo2$alphahat[1:313, "level"], col = "blue")
lines(smo2_seas[1:313] + smo2$alphahat[1:313, "level"] + smo2$alphahat[1:313, "sea_dummy1"], col = "green ") #sea_dummy1 è la stagionalità a dummy stocastiche

Trend piu stagionalità in blu. Solo trend in rosso. trend e stagionalità ogni sette giorni in verde.

mod2_ <- SSModel(c(y, rep(NA, 178)) ~ XX + css + sii +
                 SSMtrend(1, fit2$model$Q[1, 1, 1])+
                 SSMseasonal(6, fit2$model$Q[2, 2, 1], "dummy"),
                 SSMseasonal(313, fit2$model$Q[2, 2, 1], "trig", harmonics = 1:24),
                 H = exp(fit2$optim.out$par[3]))

# Poniamo anche le condizioni iniziali come nel modello precedente
mod2_$a1 <- fit2$model$a1
mod2_$P1 <- fit2$model$P1
mod2_$P1inf <- fit2$model$P1inf

# Calcoliamo lo smoothing delle variabili di stato e del "segnale"
smo2_ <- KFS(mod2_, smoothing = c("state", "signal"))

prev2 <- smo2_$muhat[test, 1]
prev2 <- if_else(prev2 < 700, prev2, prev2 + 700)


# Confrontiamo previsioni con dati reali
lin <- xts(prev2, order.by=as.POSIXct(load[test,]$Date), frequency = 6)
## Create plot of dataset
pp2 <- plot.xts(ts, screens=1,
             major.ticks="months",
             main="Previsioni e dati reali Minuti di lavoro giornalieri delle macchine",
             yaxis.right=FALSE,
             col="black")

## Add lines
pp2 <- addSeries(lin, on = 1, type = "l", col = "red", lty = 1, lwd = 1, pch = 0)

## Add legend
pp2 <- addLegend("topleft", 
              legend.names=c("real", "predicted"),
              col=c("black","red"),
              lty=c(1,1),
              lwd=c(1,1),
              ncol=2,
              bg="white")

pp2

score_UCM2_test <- sqrt( mean( (prev2-y_)^2 , na.rm = TRUE ) )
print(score_UCM2_test)

[1] 461.3332

MODELLO ARIMA CON FOURIER

library(xts)
library(forecast)
library(tseries)
library(ggplot2)
library(urca)
library(tsoutliers)
library(fpp2)
library(timeSeries)
library(KFAS)
library(tidyverse)
library(glue)
library(forcats)
library(timetk)
library(tidyquant)
library(tibbletime)
library(cowplot)
library(recipes)
library(rsample)
library(yardstick) 
library(keras)
library(forecast)
library(dplyr)

Assegnare stagionalita alla serie

ts_train <- msts(y, c(6, 26, 313.25))
ts_test <- msts(y_, c(6, 26, 313.25))

Creo la funzione per visualizzare Acf e Pacf e la applico al dataset, per trovare i picchi.

acf_pacf <- function(x, max.lag){
  par(mfrow = c(1, 2))
  Acf(x, max.lag)
  Pacf(x, max.lag)
}

Sulla base dei valori ottenuti precedentemente, sarima 1,0,1 1,1,1,6

mod1 <- Arima(ts_train, c(1,0,1), list(order = c(1,1,1), period = 6))
mod1

Series: ts_train 
ARIMA(1,0,1)(1,1,1)[6] 

Coefficients:
         ar1      ma1    sar1     sma1
      0.8609  -0.5863  0.1291  -0.9109
s.e.  0.0679   0.1001  0.0667   0.0361

sigma^2 estimated as 85223:  log likelihood=-2180.71
AIC=4371.41   AICc=4371.61   BIC=4390.05

acf_pacf(modloglik$residuals, 36)

Ora integro con i regressori di fourier:

mod3

Series: ts_train 
Regression with ARIMA(6,1,2)(1,0,1)[6] errors 

Coefficients:
          ar1      ar2      ar3      ar4      ar5     ar6      ma1      ma2     sar1     sma1
      -0.8063  -0.6723  -0.7132  -0.7374  -0.7050  0.0616  -0.5848  -0.4152  -0.3227  -0.4476
s.e.   0.1353   0.1038   0.1024   0.1050   0.1087  0.1136   0.1244   0.1242   0.1162   0.1810
          S1-6       C1-6      S2-6      C2-6     C3-6     S1-26      C1-26     S2-26    C2-26
      206.2809  -133.3336  148.8658  107.3445  99.1603  283.0346   267.8244  -20.7541  23.7697
s.e.   29.8944    29.9185   23.7046   23.6868  21.6103  737.8885  1788.5066   14.9397  15.0892
         S3-26     C3-26    S4-26     C4-26    S5-26    C5-26    S6-26    C6-26    S7-26
      -15.4524  -22.9586   3.0583  -25.7304  -2.5576  32.6571  15.0326  33.5452  -0.1282
s.e.   16.5880   16.5780  33.0044   32.6649  14.8437  14.8936  22.2930  22.3549  22.9405
        C7-26     S8-26     C8-26     S9-26    C9-26     S1-313    C1-313     S2-313   C2-313
      -8.3722  -28.8898  -17.6696  -27.9556   7.8014  -265.1517  -34.5928  -158.3469  56.5399
s.e.  22.9548   16.2744   16.3017   21.2927  21.2312     2.1281   16.1481     2.3063  16.4651
        S3-313    C3-313    S4-313   C4-313    S5-313    C5-313   S6-313   C6-313   S7-313
      -64.1502  -53.4813  -53.5860  98.3234  -37.8162  -33.2164  41.2037  25.3050  -0.1909
s.e.    2.6168   17.0226    3.0531  17.8750    3.6120   19.1108   4.3406  20.8914   5.3474
       C7-313   S8-313    C8-313    S9-313   C9-313  S10-313  C10-313  S11-313   C11-313
      20.9224  -6.5052  -84.2540  -39.3672  -9.9984  55.7634   7.6942  76.0998  -22.1824
s.e.  23.5131   6.7966   27.5662    9.1696  34.4401  13.8587  48.2444  27.6643   88.8120
        S12-313    C12-313  S13-313  C13-313  S14-313   C14-313  S15-313  C15-313   S16-313
      -229.2221  -243.7349  27.5772  -8.9120  11.1415  -43.4383  21.1590  55.4190  -23.4944
s.e.   625.2965  1827.5632  33.5548  87.5059  17.9484   40.8048  13.4487  26.1848   11.7145
      C16-313   S17-313   C17-313  S18-313  C18-313  S19-313   C19-313  S20-313  C20-313
      -5.7708  -31.9248  -34.9504  28.9982  15.9717  -4.3870  -29.2079   6.8254  29.1419
s.e.  19.4309   11.0946   15.8930  11.0540  14.0464  11.3644   13.2111  11.8802  13.0143

sigma^2 estimated as 50778:  log likelihood=-2097.59
AIC=4343.17   AICc=4390.01   BIC=4620.16

#MAPE on training
plot(ts_train, type = "l")
lines(mod3$fitted, col = "red")

legend(1, 95, legend=c("Fitted", "Train"),
       col=c("red", "black"), lty=1:1, cex=1)

score_ARIMA_train

[1] 196.9083

print(rmse1)

[1] 651.1751

print(rmse2)

[1] 856.5661

Calcolo il MAPE

PREVISIONI

autoplot(pre1)

prev_train <- forecast(pre_training, xreg=fourier(ts_train, K=c(3, 7, 20), 313), 313)
plot(prev_train)

prev_test <- forecast(pre1, xreg=fourier(ts_train, K=c(3, 7, 20), 178), 178)
plot(prev_test)

pretrain <- tail(ts(prev_train$mean+400), 174)
lower95train <- tail(prev_train$lower[,2], 174)+400
upper95train <- tail(prev_train$upper[,2], 174) +400
plot(ts(y), type = "l", main = "ARIMA train")
lines(ts(prev_train$mean), col="red")

pre <- tail(ts(prev_test$mean+400), 174)
lower95 <- tail(prev_test$lower[,2], 174)+400
upper95 <- tail(prev_test$upper[,2], 174) +400
plot(ts(y_), type = "l", main = "ARIMA test")
lines(ts(prev_test$mean), col="red")

#export data
write.csv(pre, file=paste(folder_risultati, "/prediction_fourier.csv", sep=""), sep = ";", dec = ",", row.names=FALSE)

attempt to set 'sep' ignoredattempt to set 'dec' ignored

write.csv(lower95, file=paste(folder_risultati, "/lower_fourier.csv", sep=""), sep = ";", dec = ",", row.names=FALSE)

attempt to set 'sep' ignoredattempt to set 'dec' ignored

write.csv(upper95, file=paste(folder_risultati, "/upper_fourier.csv", sep=""), sep = ";", dec = ",", row.names=FALSE)

attempt to set 'sep' ignoredattempt to set 'dec' ignored